%% master_counterfactuals_CES_idiosyncraticprefshock.m

clear

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author: Produced for JdG,AL,IM
%
% Date: July 2019
%
% Purpose: This file runs the counterfactual exercises after the
% main_baseline.m has been run. Calibrate idiosyncratic productivity shocks
% that are equivalent in size to firms' exposure to the foreign shock
% 
% Based on one file: Step_4_Algorithm_counterfactual.m, which uses new
% baseline input data, MAT/Step_4_setup_`alphaT''rhoT'_`shockT'.mat
% as initial level variables, and then shocks the system and computes hat
% variables for wages, prices, exports, etc.
%
% NEW: for CES setup. Have streamline file that has to be called up and
%                     various variables that have to be used.
% Output: MAT/Step_4_results_`alphaT''rhoT'`lambdaT'`etaT'_`shockT'.mat,
%         MAT/Step_5_results_`alphaT''rhoT'`lambdaT'`etaT'_`shockT'.mat
%         MAT/Elasticity_output_`alphaT''rhoT'`lambdaT'`etaT'_`shockT'.mat
%
% Last Updated: 07/09/2019 by JDG/IM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global alphaT rhoT lambdaT etaT shockT rho sigma sigmaT lambda eta ...
    tol maxit Pi_hat_n D_hat_n Xi_hat_mnjf Xi_hat_mnj a_hat_f a_hat ...
    varphi sL_n sPi_n sD_n N J FI_J france oneminus_alpha_vector ...
    labour_share_PWT countrysecD firmsecD_sorted start start_sorted ...
    sum_by_country_dummy ISO alpha_WIOD_j alpha_WIOD_francej sector_algorithm 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% alphaT: 
%
% 1 = PWT labor share for all countries and Labor share for French firms 
%     are at the COUNTRY level.
%
% 2 = PWT labor share for all countries and WIOD scaled labor share for 
%     French firms. Labor share for French firms are at the FIRM level.
%
% 3 = WIOD labor share for all SECTORS with the labor share at the SECTOR 
%     level for French firms.
%
% 4 = WIOD labor share for all SECTORS with the labor share at the FIRM 
%     level for French firms 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Loading in the previous data to get some useful variables

alphaT = 4;
rhoT = '3';
lambdaT = '1.000001';
etaT = '1.000001';
lambda = 1.000001;
eta = 1.000001;
sigmaT='1.5';
%rescale_factor=0.039023637/0.178732791359661;


fname = strcat('MAT/Step_2_firm_data_alpha',num2str(alphaT),'.mat');
load(fname,'FI_J', 'firmsecD_sorted', 'start_sorted');
load('MAT/WIOD.mat','pi_c_mnj','countrysecD','france','J','N','start','sum_by_country_dummy','gamma','alpha_WIOD_nj');


rho = 3*ones(J,1);
sigma = 1.5*ones(J,1);

%% SET PARAMETERS VALUES

% Frisch elasticity for GHH preferences

varphi=3;

% share of labour in total final consumption (need to make assumption or collect data)

sL_n = ones(1,N)*0.8; 

% share of profits in total final consumption (should have this data)

sPi_n= ones(1,N)*0.15; 

% share of trade deficit in total final consumption (should have this data)

sD_n= ones(1,N)*0.05; 

%% Baseline: no shocks

% Pi_hat: Change in profits per country
    
Pi_hat_n=ones(1,N);
    
% D_hat: Change in trade deficit
    
D_hat_n = ones(1,N);
    
% Taste shock - for France only
    
Xi_hat_mnjf=ones(FI_J,N);
    
%Taste shock - for sectors only (countries except France)
    
Xi_hat_mnj= ones(N*J,N);
    
% Productivity shock, a is firm specific - France only
    
a_hat_f= ones(FI_J,1);
    
% Productivity shock, a is for sectors - countries except France
    
a_hat = ones(J*N,1);

%% creates an index that captures each sector-decile set of firms
eval([strcat('load MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),'_sigma',sigmaT,'/Step_5_output_alpha',num2str(alphaT),...
    '_rho',rhoT,'_WORLDprod_p10','_approx.mat')]);

ps=1;
PS=[];
for j=1:J
    perc_sec = zeros(size(VA_fj_0(start_sorted(j):start_sorted(j+1)-1,1)));
    p = prctile(full(VA_fj_0(start_sorted(j):start_sorted(j+1)-1,1)),99);
    perc_sec(VA_fj_0(start_sorted(j):start_sorted(j+1)-1,1)>=p & perc_sec ==0) = 1;
    p = prctile(full(VA_fj_0(start_sorted(j):start_sorted(j+1)-1,1)),90);
    perc_sec(VA_fj_0(start_sorted(j):start_sorted(j+1)-1,1)>=p & perc_sec ==0) = 9;
    perc_sec(perc_sec ==0) = 90;
    PS = [PS;perc_sec];
end

shockT = 'DOMidioprod'
             
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Homogeneity function
% %
% % We need to homogeneize (1-alpha)gamma given the structure of the model,
% % along with the labor and intermediate shares in the data, for French
% % firms. Having done this, we can then run the model as for the
% % heterogeneous case.
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% % X_nkif: Multiplies the expenditure matrix for france with the firm dummy
% %         to create a [F x N matrix]
% 
% X_nkif = firmsecD_sorted*X_mnj(start(france):start(france+1)-1,:);
% 
% % pi_X_nkif: Multiplying firm shares of expenditures (pi_nkj) by the
% %            expenditure matrix that has been scaled up by the firm 
% %            sector dummy
% 
% pi_X_nkif=pi_mnjf.*X_nkif;
% 
% % pi_X_nif: This is the sum over country k, should result in a [Fx1] vector
% 
% pi_X_nif=sum(pi_X_nkif,2);
% 
% %one_minus_pi_l_f: Labor shares for each firm = 1-pi_l_f
% 
% one_minus_pi_l_f=1-pi_l_f;
% 
% % Multiplying the previous summed equation by this labor share. This is
% % done before multiplying by pi_M_f, which will come later.
% 
% pi_X_laborshare_nif=pi_X_nif.*one_minus_pi_l_f;
% 
% clear pi_X_nkif  X_nkif
% 
% pi_X_laborshare_nif_repmat = repmat(pi_X_laborshare_nif, [1 J*N]);     
% 
% gamma_pi_X_laborshare_nif = pi_X_laborshare_nif_repmat.*pi_M_f;
% 
% % Sum over firms. Should give a J*N matrix, which is numerator of
% % homogeneous (1-alpha)*gamma
% 
% numerator = firmsecD_sorted'*gamma_pi_X_laborshare_nif;
% 
% % Sum exports over destination country per sector for france for numerator
% % 32*1 matrix, which must be expanded per destination market
% 
% X_fra = X_mnj(start(france):start(france+1)-1,:);
% den = sum(X_fra,2);   
% denominator = repmat(den,[1 J*N]);
% 
% clear pi_X_laborshare_nif_repmat gamma_pi_X_laborshare_nif X_fra 
% 
% % Create homogeneous (1-alpha)gamma
% 
% one_minus_alp_gam_hom = numerator./denominator;
% 
% clear numerator denominator;
% 
% % Create homogeneous alpha
% pi_X_laborshare_nif=pi_X_nif.*pi_l_f;
% numerator = firmsecD_sorted'*pi_X_laborshare_nif;
% pi_l_hom = numerator./den;
% 
% clear numerator den pi_X_laborshare_nif;
% % Create homogeneous gamma given that alpha_hom = alpha_baseline at
% % country-sector level
% 
% pi_l_hom_check = pi_l(start(france):start(france+1)-1,:);
% corr(pi_l_hom,pi_l_hom_check) 
% pi_M_hom = one_minus_alp_gam_hom./(1-repmat(pi_l_hom,[1 N*J]));
% 
% % Now fill in firm levels given homogeneous values and sectors firms in
% 
% pi_l_f_hom = firmsecD_sorted*pi_l_hom;
% pi_M_f_hom = firmsecD_sorted*pi_M_hom;
% 
% clear denominator numerator one_minus_alp_gam_hom pi_l_hom pi_M_hom;
%          
% 
% pi_M_f_foreign = pi_M_f_hom;
% pi_M_f_foreign_g_P = pi_M_f_foreign.* repmat(transpose(g_P),[FI_J,1]);
%         %pi_M_f_foreign(:,start(france):start(france+1)-1)=0; 
% sum_pi_M_f_foreign_g_P = sum(pi_M_f_foreign_g_P,2);
% one_minus_alpha = 1-pi_l_f_hom;
% one_minus_alpha_sum_gamma_foreign_g_P = one_minus_alpha .*sum_pi_M_f_foreign_g_P;
% exposure = one_minus_alpha_sum_gamma_foreign_g_P;
% exposure_deciles = exposure.*Quantiles_Dummy; 
% clear P_hat;


for sim=[1 9 90]

    Pi_hat_n=ones(1,N);
    D_hat_n = ones(1,N);
    Xi_hat_mnjf=ones(FI_J,N);
    Xi_hat_mnj= ones(N*J,N);
    
    a_hat = ones(J*N,1);
    a_hat_f=ones(FI_J,1);
    a_hat_f(PS==sim,1)=.9091;
    
        %% Running the Algorithm

        disp('BEGINNING COUNTERFACTUAL EXERCISE')

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
        % Step_4_Algorithm: same algorithm as Step_4_Algorithm_baseline.m, but now
        % we are able to feed in the model the calculated t+1 varialbes as t,
        % starting at the new equilibrium with no wedges
    load MAT/WIOD.mat
    eval([strcat('load MAT/Step_2_firm_data_alpha',num2str(alphaT),'.mat')]);
    eval([strcat('load MAT/rho',num2str(rhoT),'eta',num2str(etaT),'lambda',num2str(lambdaT),'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_3_alpha',num2str(alphaT),...
    '_rho',rhoT,'.mat')]);
    eval([strcat('load MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_4_setup_alpha',num2str(alphaT),...
    '_rho',rhoT,'_approx.mat')]);


         Step_4_Algorithm_counterfactual_CES_Hom_alphagamma_approx
% 
%         % Save the counterfactual output so it can be used for post analysis
% 
         outputtable = [' X_mnj X_hat_mnj pi_mnjf1 pi_mnj1 ' ...
                  'pi_c_mnj1 pi_c_nj1 pi_M1 pi_l1 pi_M_f1 pi_l_f1 '...
                  'PC_hat_n1 PC_n P_hat_n P_hat_nj P_hat_mnj w_hat ' ...
                  'a_hat a_hat_f Xi_hat_mnj Xi_hat_mnjf '...
                  'final_goods1 intermediates1'] ;
%         %  
         eval([strcat('save MAT/rho',num2str(rhoT),'eta',num2str(etaT),'lambda',num2str(lambdaT),'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_4_results_alpha',num2str(alphaT),...
                  '_rho',rhoT,'_',shockT,'_approx_Hom_Percentile',num2str(sim),'.mat',outputtable,' -v7.3;')]);   
%         %     
         display('STEP 4 COUNTERFACTUAL COMPLETED SUCCESFULLY')

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


            % Step_5_counterfactual_output: Allows the user to easily analyse the
            %      counterfactual exercises produced before. Runs the CounterOutput.m
            %      function after specifying alpha type (labor share), gamma
            %      calibration and the shock scheme to be analysed.

            %% Load data

            
            % Time 0 data

            eval([strcat('load MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,...
    '_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_4_setup_alpha',num2str(alphaT),...
     '_rho',rhoT,'_approx.mat')]);

            % Time 1 data

            eval([strcat('load MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,...
     '_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_4_results_alpha',num2str(alphaT),...
     '_rho',rhoT,'_',shockT,'_approx_Hom_Percentile',num2str(sim),'.mat')]);
 
        Step_5_counterfactual_output_CES_Hom_alphagamma_approx

            % Save the counterfactual output so it can be used for post analysis
        outputtable = [' VA_hat_n VA_hat_mj VA_hat_mnj VA_hat_fj VA_hat_fnj '...
          'RGDP_def_hat_n RGDP_def_hat_fj RGDP_cpi_hat_n RGDP_cpi_hat_fj '...
          'RGDP_doubledef_hat_n RGDP_doubledef_hat_fj '...
          'real_wage_cpi_hat real_wage_def_hat real_wage_doubledef_hat '...
          'real_wage_income_cpi_hat_n real_wage_income_def_hat_n real_wage_income_doubledef_hat_n '...
          'TFP_hat_def_n TFP_hat_cpi_n TFP_hat_doubledef_n TFP_hat_cpi_france_j '...
          'imports_over_GDP_all imports_over_GDP_france_sector '...
          'exports_over_GDP_all exports_over_GDP_france_sector '...
          'GDP_deflator_n DoubleDeflation_deflator_n P_hat_n VA_fj_0 VA_mj_0'] ;
        %  
        sigmaT = '1.5';
             eval([strcat('save MAT/rho',num2str(rhoT),'eta',num2str(etaT),'lambda',num2str(lambdaT),'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_5_output_alpha',num2str(alphaT),...
                 '_rho',rhoT,'_',shockT,'_approx_Hom_Percentile',num2str(sim),'.mat',outputtable,' -v7.3;')]); 

        %     
             display('STEP 5 COMPLETED SUCCESFULLY')

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            % Step_6_Elastisticity_output: perform covariance decomposition at firm and
            %        sector level

            % Load files: only needed if run separately from Step 5

            eval([strcat('load MAT/rho',num2str(rhoT),'eta',num2str(etaT),'lambda',num2str(lambdaT),'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_5_output_alpha',num2str(alphaT),...
                 '_rho',rhoT,'_',shockT,'_approx_Hom_Percentile',num2str(sim),'.mat;')]); 

            % Shock size

          delta = 0.1;

           Step_6_Elasticity_output_CES

       outputtable = [' e_Y_FRf_cpi m_e_f_cpi cov_e_f_cpi check_f_cpi '...
           'e_Y_FRs_cpi m_e_mjFRs_cpi cov_e_mjFRs_cpi check_mjFRs_cpi '...
           'e_Y_s_cpi m_e_mj_cpi cov_e_mj_cpi check_mj_cpi '...
           'e_Y_FRf_def m_e_f_def cov_e_f_def check_f_def '...
           'e_Y_FRs_def m_e_mjFRs_def cov_e_mjFRs_def check_mjFRs_def '...
           'e_Y_s_def m_e_mj_def cov_e_mj_def check_mj_def '...
           'e_Y_FRf_doubledef m_e_f_doubledef cov_e_f_doubledef check_f_doubledef '...
           'e_Y_FRs_doubledef m_e_mjFRs_doubledef cov_e_mjFRs_doubledef check_mjFRs_doubledef '...
           'e_Y_s_doubledef m_e_mj_doubledef cov_e_mj_doubledef check_mj_doubledef'] ;

       eval([strcat('save MAT/rho',num2str(rhoT),'eta',num2str(etaT),'lambda',num2str(lambdaT),'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_6_output_alpha',num2str(alphaT),...
                 '_rho',rhoT,'_',shockT,'_approx_Hom_Percentile',num2str(sim),'.mat',outputtable,' -v7.3;')]); 

             display('STEP 6 COMPLETED SUCCESFULLY')

        disp('COUNTERFACTUAL EXERCISE COMPLETED SUCCESFULLY')
        

end
